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The aim of this work is the description of the chain formation phenomena observed in colloidal 
suspensions of superparamagnetic nanoparticles under high magnetic fields. We propose a new 
methodology based on an on-the-fly Coarse-Grain (CG) model. Within this approach, the coarse 
grain objects of the simulation and their dynamic behavior are not fixed a priori at the beginning 
of the simulation but rather redefined on-the-fly. The motion of the CG objects (single particles or 
aggregates) is described by an anisotropic diffusion model and the magnetic dipole-dipole interac- 
tion is replaced by an effective short range interaction between CG objects. The new methodology 
correctly reproduces previous results from detailed Langevin Dynamics simulations of dispersions 
of superparamagnetic colloids under strong fields whilst requiring an amount of CPU time orders 
of magnitude smaller. This substantial improvement in the computational requirements allows the 
simulation of problems in which the relevant phenomena extends to time scales inaccessible with 
previous simulation techniques. A relevant example is the waiting time dependence of the relax- 
ation time T2 of water protons observed in magnetic resonance experiments containing dispersions 
of superparamagnetic colloids, which is correctly predicted by our simulations. Future applica- 
tions may include other popular real-world applications of superparamagnetic colloids such as the 
magnetophoretic separation processes. 

PACS numbers: 05.10.-a, 82.70.Dd, 87.15.nr, 47.65.Cb 



I. INTRODUCTION 

In recent years, work in coarse grain models for the 
description of soft matter and biomolecular systems is 
experiencing a remarkable outburst [T]. The reason is 
that the description of these systems at experimentally 
relevant time and length scales requires inclusion of phe- 
nomena occurring at very different scales. The objective 
of coarse grain (CG) models is thus to retain sufficient 
molecular or nanoscale detail and yet remain amenable 
of simulation up to macroscopic time scales. Many ap- 
proaches have been developed to construct CG models of 
different kinds of soft matter systems. For example, in 
the case of polymers, there is a long tradition of using CG 
models and the field is sufficiently mature so that there 
are systematic and rigorous approaches to build up CG 
models from accurate atomistic descriptions [5] . Also, in 
the field of biomolecular simulations, there are important 
developments such as the MARTINI force field |3_ which 
allow the simulation of difficult problems such as the be- 
havior of lipid vesicles or protein folding at millisecond 
or even larger time scales. New advances include also 
simulation packages specially designed for CG models of 
soft matter such as ESPResSo [I]. 

Our interest here is the development of an improved 
CG model for a specific problem which is still difficult 
to simulate, namely the assembly of superparamagnetic 
colloids under strong magnetic fields. Superparamag- 
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netic colloids are made of small nanoparticles of magnetic 
material (typically 5-10 nm iron oxide nanocrystals) em- 
bedded in a nonmagnetic matrix (typically polymers or 
silica) |5_ . These particles have no magnetic dipole in ab- 
sence of magnetic field but they develop very high mag- 
netizations in the presence of a magnetic field, similar to 
those obtained with ferromagnetic materials. This highly 
tunable response and the possibility to functionalize their 
surface make these materials very interesting for appli- 
cations such as capture and removal of biomolecules and 
pollutants, NMR contrast agents, and many others (SHE]- 

Our work is motivated by the difficulties encountered 
in modeling different kinds of real experimental situations 
involving superparamagnetic colloids. A relevant exam- 
ple is provided by the experiments by Chen et al. [5] of 
a dispersion containing superparamagnetic colloids de- 
signed as contrast agents for magnetic resonance imag- 
ing (MRI). In these experiments, a strong uniform mag- 
netic field was applied to the dispersion. It was observed 
that the transverse relaxation time T 2 of protons in wa- 
ter changed with time, an effect which was attributed to 
the formation of chains of superparamagnetic colloids. In 
fact, the kinetics of chain formation was estimated from 
these experiments, spanning time scales from 10 to 10 3 
seconds or more. Another relevant example is magne- 
tophoresis |10| . which is the motion of magnetic parti- 
cles under an external magnetic gradient. Experimen- 
tal evidence shows that the formation of chains induced 
by the external field speeds up the magnetophoresis pro- 
cess |11J . which is orders of magnitude faster than that 
observed in absence of chain formation [12]. It is worth 
noting that in these experiments, chains dissolve almost 
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immediately after removal of the external field, as should 
be expected since superparamagnetic nanoparticles have 
no dipole in absence of magnetic field. In this respect, 
these systems are very different from the widely studied 
(and simulated) dispersions of dipolar particles, which 
are able to form structures in the absence of an external 
magnetic field due to the interaction of their permanent 
dipoles |13fTT5] , 

The standard approach for simulation of chaining pro- 
cesses in magnetic colloids is the use of Langevin Dy- 
namics (LD) simulations (see for example [ 1 6] [TS] ) . This 
technique allows the inclusion of particle-particle inter- 
actions, thermal noise and the friction due to the fluid. 
The resolution of the simulation technique is typically 
in the nanoseconds scale. Simulation runs up to a 
few seconds are possible, but they are highly intensive, 
requiring the use of parallel computing during several 
weeks [18 . These CPU requirements make this simula- 
tion technique unsuitable for the study of the problems 
mentioned above. 

The need to account for microscopic time and length 
scales but also reach macroscopic time scales at low com- 
putational cost has motivated us the development of a 
new simulation strategy based on an on-the-fly CG pro- 
cedure. The methodology which will be developed in 
this paper is a generalization of two procedures pro- 
posed in previous works: the method proposed by Miguel 
and Satorras jT3] to study aggregation processes and the 
method proposed by Schaller et al. [20] to study magne- 
tophoresis. 

In the methodology proposed here, one starts by simu- 
lating the motion and interaction between individual col- 
loids. As the simulation advances, colloids form chains 
due to the magnetic dipole-dipole attraction induced by 
the high external magnetic field. The motion of each 
particle inside a chain is not simulated explicitly. In 
our methodology, these chains are considered individ- 
ual coarse grain (CG) objects which move following cer- 
tain effective rules and interact (and possibly aggregate) 
with other CG objects or single individual particles. In 
this way, the CG objects of the simulation are not fixed 
a priori at the beginning of the simulation but rather 
redefined on-the-fly. Thus, we adjust the resolution of 
the calculations during the simulation run, allowing for 
the possibility of much longer simulation runs requir- 
ing less computer power. Preliminary simulation results 
and comparison with experiments, presented in a previ- 
ous work [5], demonstrated the feasibility and utility of 
our novel approach. Here we will discuss in detail the 
physical basis of the model, the simulation methodology 
and detailed comparison with more standard Langevin 
simulation techniques. All simulations of our model were 
performed employing the MagChain program, a CH — h ap- 
plication developed in-house, which is freely available for 
use of researchers. The code, its documentation and us- 
age examples can be found available for download at our 
web page [21J. 

The paper is organized as follows. In Section II, we 



describe the modeling of the system under study and the 
simulation technique. In Section III we validate the new 
methodology by (i) comparing our results with those ob- 
tained employing standard LD simulations (ii) discussing 
the effect of choosing other approximations for the dif- 
fusion coefficients and the effective interaction of the 
CG objects and (iii) discussing the applicability of our 
methodology comparing with experimental results. The 
conclusions are presented in section IV and some techni- 
cal issues are detailed in the Appendix. 



II. FORMULATION OF THE MODEL 

The system which we are interested to describe is a 
colloidal dispersion of N p superparamagnetic spherical 
particles of diameter d in a volume V and volume fraction 
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In absence of external magnetic field, the particles have 
no magnetic dipole and there is no formation of chains 
(no aggregation induced by the magnetic field). In the 
presence of a magnetic field H, the superparamagnetic 
particles acquire a certain magnetization M(H). Since 
we are particularly interested in the case of very strong 
magnetic fields (as in the experiments of Ref. [5J for ex- 
ample), we consider that the particles have a magnetic 
dipole moment m s (corresponding to the saturation mag- 
netization M s of the particles) pointing in the direction 
of the applied magnetic field (which we will take as the z 
axis). The strength of the magnetic interaction between 
particles as compared with thermal energy can be char- 
acterized by the magnetic coupling parameter T defined 



Mo 



(2) 



The behavior of superparamagnetic colloids under exter- 
nal fields is controlled by the values of these two pa- 
rameters (0o and r). In this paper, we are interested 
in situations (values of <po & n d T) in which the external 
field induces formation of linear chains of colloids which 
grow irreversibly with time. Irreversible growth of linear 
chains has been found in simulations and experiments 
investigating the ranges of T between 40 and 3xl0 3 and 
4>o < 0.15 |16tfTH] . However, other structures are found 
at different ranges of <po and T. For lower values of T, 
an equilibrium state is possible, in which colloids aggre- 
gate in linear (non branched) chains with an equilibrium 
length given by \/ 0oe r_1 [IB]. In the opposite situation 
of larger values of </>o and T, different aggregate structures 
can be found, including thick chains (obtained from lat- 
eral aggregation of linear chains) , bundles and more com- 
plex fibrous structures [17] . All these more complex situ- 
ations, different from irreversible growth of linear chains, 
will be left for future extensions of the model. 
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As key ingredients to retain the underlying physics of 
irreversible chain growth, we consider both the diffusive 
motion of particles and chains and their respective mag- 
netic and steric interactions. The main approximations 
of our model will be to ignore the details of the particles 
forming the chains and to replace the actual magnetic 
dipole-dipole interaction by an effective short-range in- 
teraction, less demanding from the computational point 
of view. 

Our model to study the kinetics of chain formation in 
these systems consists of CG objects which are chains 
made of s particles, including the case s = 1 which cor- 
responds to a single particle. The first ingredient of the 
model is the diffusion coefficient of the CG objects. For 
single, isolated particles (s = 1) we have: 



k B T 
3Tri]d ' 



(3) 



where rj is the viscosity of the fluid. A chain containing 
s > 1 particles exhibits anisotropic diffusion, character- 
ized by a diffusion coefficient Ds in the direction parallel 
to the long axis of the chain and Djr in the directions 
perpendicular to the long axis. There are several pos- 
sibilities for the analytical form of these diffusion coef- 
ficients, depending on the exact geometry assumed for 
the chains and the degree of approximation of the cal- 
culation. Here, in order to keep the model as simple as 
possible, we consider the following expressions valid for 
elongated objects (slender body theory [22J): 



Di 2s 



^- = -[ln(2s) 
Di 4s L V ' 
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Strictly speaking, Eqs. Q and J5| are valid only for large 
s. Therefore, we employ Eqs. Q^and ^ for chains with 
s > 2 and use a simple interpolation between the diffu- 
sion coefficients corresponding to CG objects with s = I 
(Eq.[3) and s — 3 (Eqs. Q and (JS)) for chains with s = 2. 
Such a choice gives results indistinguishable from those 
provided by more sophisticated and accurate expressions 
of the diffusion coefficients (see Section IIIB). 

The second ingredient of the model is the definition of 
the effective interaction between CG objects. A CG ob- 
ject of length s interacts with other CG objects through 
an excluded volume interaction (hard core) correspond- 
ing to a cylinder of length s x d and diameter d. They 
also interact through dipole-dipole interactions. In or- 
der to simplify and speed up the simulations, we have 
replaced the actual dipole-dipole magnetic interaction 
between colloids by an effective, short range interaction 
between the CG objects. This interaction is defined as 
follows. For a given CG object, we define two spheri- 
cal attractive regions of radius r a (s) (which depend on 
the length of the chain s) located at the two ends of the 



chain. As illustrated in Fig.[T] these regions are designed 
to mimic the region at which the magnetic attraction be- 
tween a chain of particles (magnetized in the z direction) 
and an incoming test dipolar particle is equal or stronger 
than the thermal energy k B T. The values of r (s) are 
calculated by finding the distance in the z axis at which 
the magnetic interaction energy E mag between a chain of 
s particles and a single test particle is equal to —k B T. 
Therefore, r a (s) is given by the solution of: 



E 



mag 



k B T 
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n ^(2r a (s)/d+l/2 + n)s 
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The results of Eq. (|6| for different values of T are also 
shown in Fig. [T] Once we have defined the range of the 
interaction, we need to define the strength of this inter- 
action. In order to keep our model as simple as pos- 
sible, we simply assume that all events in which a CG 
object enters into the interaction region of another CG 
object will lead to instantaneous aggregation. This rule 
has been employed previously in the interpretation of ex- 
perimental results and it has been suggested by direct 
observations of chain formation under a microscope (see 
Refs. |23l I24j ). As we will show in the next section, this 
rule reproduces correctly the results obtained from LD 
simulations in which the magnetic interaction is com- 
puted accurately. The sensitivity of the results to the 
choice of r a will be also discussed in Section IIIB. 

Once the basic ingredients for the model (rules for mo- 
tion and interaction) are defined, it is necessary to specify 
the algorithm for the numerical solution of the model. In 
our case, the diffusive motion of the CG objects is sim- 
ulated using the brownian dynamics technique [25]. At 
each time step At a random displacement in each direc- 
tion is generated with a gaussian distribution with zero 
mean and variance 2D s At, where D s is the diffusion co- 
efficient of the object (single particle or chain) in the 
direction of motion (x, y or z). Also, at each time step 
the distances between CG objects are checked in order 
to detect penetration of a CG object inside the region of 
aggregation of another CG object, as explained above, or 
to detect possible overlaps between them. In the case of 
aggregation of two CG objects, a new CG object is cre- 
ated (and the two previous CG objects are erased from 
the simulation) with length s, obtained from adding the 
lengths of the two aggregating chains and located at the 
center of mass of the aggregating CG objects. In the 
case of overlap between two CG objects without pene- 
tration into the aggregation region, we consider that the 
two chains collide. In this case, the moving CG object 
is placed in contact with the other one (without overlap- 
ping) at the collision coordinates defined by the trajec- 
tory previously followed (see Fig. [2]) . Finally, it should 
be noted here that the selection of an appropriate time 
step At for the simulation is a crucial issue. A detailed 
discussion on the selection of At is given in the Appendix. 

Hence, a typical simulation run is as follows. The sim- 
ulation starts from a pre-equilibrated system containing 
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FIG. 1: (Color online) (a) 2D-map corresponding to the interaction energy between an incoming test dipolar particle and 
a chain-like aggregate formed by 5 colloids with a magnetic coupling of F=40. The black dashed line delimites the region 
with (E < —hsT). The region excluded by the finite size of the 5 spheres is shown in blue (the interaction energy map is 
not evaluated inside this region), (b) Sketch (top view) of the attraction model implemented in the MagChain code. Each 
CG object has two attraction zones, modeled as a sphere of radius r a (Eq. (JgJ) tangent to the edge of the aggregate. Any 
particle entering into these zones will immediately aggregate forming a longer chain, (c) Dependence of the radius (r a ) of the 
attraction regions on the aggregate size for two different values of the coupling parameter, T=40 (open symbols) and T=247 
(solid symbols). The attraction radius increases abruptly for short chains and tends to a constant value for longer chains. All 
the distances are expressed in terms of the diameter d of the colloid. 
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tail in the following section. As a simulation output, we 
obtain the number of chains containing s colloidal par- 
ticles at time t, n s {t). During the simulation, we also 
monitor the time evolution of the average number of col- 
loidal particles in a chain (N)(t) defined as in |161 119j : 



(N(t)) 



(7) 



and the probability of finding an aggregate of size s at a 
given time, defined as: 



p(s;t) 



njf) 



(8) 



FIG. 2: (Color online) Sketch corresponding to the scheme 
applied to avoid the the overlap between CG objects during 
simulations. If the random displacement performed on object 
B produces an overlap between B and another CG object (A), 
the moving aggregate (B) is placed in contact with the second 
aggregate (A) according to the trajectory followed during the 
random displacement. 



N p colloids (CG objects with s = 1). As the simulation 
goes on, colloids aggregate and chains with increasing val- 
ues of s appear. Consequently, the number of CG objects 
of the simulation decreases with time and the simulation 
speeds up as the time advances, as we will discuss in de- 



III. VALIDATION AND APPLICATION OF THE 
MODEL 

A. Comparison with Langevin Dynamics 
simulations 

Our objective in this section is to compare the perfor- 
mance and results obtained using the model described 
in Section II with standard LD simulations of the same 
system. Briefly stated, Langevin Dynamics simulations 
consist on solving the Newton equations of motion for 
each particle taking into account external forces, the in- 
teraction forces between particles (magnetic and steric), 
the viscous drag from the solvent and a stochastic force 
arising from the thermal noise due to the fact that the 
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TABLE I: Characteristics of the colloidal dispersions of superparamagnetic particles simulated with CG and LD techniques. 
(fio and r are defined by Eqs. ([I]) and pj, p p is the density of a single colloid, d is its diameter and D\ its diffusion coefficient. 
T is the temperature of the dispersion and 77 the viscosity of the solvent (water) at this temperature. 
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00 


Pp [g/cm 3 ] 


d [nm] 


T [K] 


r\ [Pa ■ s[ 


Lh. [m 2 /s] 


Case 1 


40 


5.23xf0" 4 


f.O 


100 


300 


l.OxlO" 3 


4.39 xl0~ 12 


Case 2 


247 


4.64xfCT 6 


3.1 


88 


310 


0.692xl0~ 3 


7.46 x 10~ 12 



TABLE II: Set of parameters used for numerical integration in the Coarse Grain (CG40 and CG247) and the Langevin Dynamics 
(LD40 and LD247) simulations. N p is the number of particles in the simulation, L z and L x — L y are the sizes of the simulation 
box (in units of particle diameters) in the directions parallel and perpendicular to the magnetic field, respectively (periodic 
boundary conditions were employed in all simulations). At is the time step and tf is the total simulated time. We also indicate 
the total amount of CPU time employed in the calculation, calculated as the number of cores used times the total elapsed time. 
In all our calculations we have used a 8-Core AMD Opteron Magny Cours 6136 processor. 
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CPU cost 


LD40 


Case 1 


8000 


200 


200 


1.02xl0" 9 


2.04 


8 


998h 


LD247 


Case 2 


4000 


767 


767 


3.06xl0 -9 


6.12 


8 


866h 


CG40 


Case 1 


8000 


512 


128 


2.280 x 10~ 4 


5 


1 


25h 


CG247 


Case 2 


8000 


1534 


767 


1.038 x 10" 4 


1000 


1 


24h 




FIG. 3: (Color online) Snapshots from simulations with F = 40 (Case 1 of Table ffl). (a) Initial configuration of a simulation 
(t = 0) (b) Snapshot from Langevin Dynamics simulations (LD40) at t — 0.28s. Note that the simulation resolves the individual 
particles building up the chains, (c) Snapshot from Coarse Grain simulations at t — 0.28s. Now the chains are the CG objects, 
individual particles are no longer considered once they form part of a chain. Chains are colored according to their length for an 
easier visualization. Left and center images were created using VMD [26J. Right image was created using our own visualization 
software available in the web |21| . 



system is at a given temperature T. This comparison be- 
tween our simplified CG method and more detailed LD 
simulations will help to clarify the validity of the approx- 
imations introduced in our model, as described in the 
previous section. In order to perform a significant com- 
parison between the new procedure and the standard LD 
simulation technique, we have selected two cases with 
very different magnetic coupling T which were studied in 
previous works. The details for these systems are sum- 



marized in Table |T] and were denoted as Case 1 and Case 
2. 

Let us consider first Case 1, which corresponds to a dis- 
persion of 100 nm superparamagnetic colloids at a vol- 
ume fraction </>o = 5.23 x 10~ 4 which have a magnetic 
coupling parameter T = 40 at saturation (i.e. under 
strong magnetic fields). This system was studied em- 
ploying LD simulations in Ref. |18| by using the standard 
LAMMPS simulation package [57] (version 21May2008). 



10" 



z 10 1 
V 



CG40 —€>•-- 
LD40 — 
CG247 — e— 
LD247 



<tt I 



10" fe: 
1(T 




L0"- 



10" 



10"' 10" 10' 
time (s) 



L0 3 



FIG. 4: (Color online) Time evolution of the average num- 
ber of particles (N(t)) in a chain, Eq. Q. Comparison be- 
tween the results obtained from Langevin Dynamics (solid 
symbols) and Coarse Grain (open symbols) simulations for 
the two different systems studied. Circles correspond to 
Case 1 (r=40,<?(>o=5.23xl(n 4 ) and squares to Case 2 (r=247, 
<?!>o=4.64xl0- 6 ). 




FIG. 5: (Color online) Comparison between the probability 
distribution to find an aggregate of size s (defined in Eq. (J1J) 
at t=ls obtained from LD40 and CG40 simulations. 



Now, we will compare these published results obtained 
with the standard LD technique with our CG methodol- 
ogy described in the previous section. These simulations 
will be denoted as LD40 (the Langevin Dynamics case) 
and CG40 (our Coarse Grain methodology) . The param- 
eters employed in the numerical algorithm are given in 
Table [TT] For completeness, we also give the parameters 
employed in the LD simulations. It should be noted that 
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FIG. 6: (Color online) Comparison between the probability 
distribution to find an aggregate of size s at t=5s obtained 
from LD247 and from CG247. 



the LD simulations require a small time step (of the or- 
der of the ns). This small time step is needed in order to 
avoid a simulation crash in simulations involving chains 
of colloids, since the motion of colloids inside a chain in- 
volves very small displacements which need to be resolved 
with high precision. In our CG methodology (which do 
not consider the structure of chains), we can use much 
larger time steps as shown in Table [TTJ A detailed discus- 
sion on the selection of the time step in our methodology 
is given in the Appendix. 

Typical snapshots of the simulation are shown in Fig. [3] 
and the results for (N(t)) are shown in Fig. [4] The 
snapshots illustrate the different resolution employed in 
the CG40 and LD40 simulations. As seen in the snap- 
shots, the LD40 simulation resolves the individual parti- 
cles making up the chains whereas the chains are struc- 
tureless in the CG40 simulation. It should be noted 
that the chains obtained in the LD40 simulation are 
almost perfectly linear and are not significantly differ- 
ent from the coarse-grain objects of the CG40 simula- 
tion. As shown in FigEJ the values of (N(t)) obtained 
from both simulations (CG40 and LD40) are in excellent 
agreement. For example, at t = Is, the mean aggregate 
size for the CG40 simulation is (iV)c6?=12.10 and the 
value calculated from LD simulations is almost identical, 
(N) 1,0=12.14. Therefore, we can conclude that the sim- 
plifying approximations included in the new methodology 
(particularly those regarding the calculations of particle- 
particle magnetic interaction) do not affect the average 
size of chains. 

We have performed a more detailed comparison be- 
tween both approaches by comparing the distribution 
of chains of size s at certain times. In Fig. [5] we com- 
pare the corresponding probability distribution (defined 
in Eq. ([§])) at t =ls obtained from LD40 and CG40 sim- 
ulations. The agreement between both results after Is 
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is remarkable, and only slight differences are observed. 
As shown in Fig. [5j the distribution of chain sizes is very 
broad, with significant probabilities of finding chains well 
above and well below the average length (including iso- 
lated particles). 

Now, let us consider the case denoted as Case 2 in Ta- 
ble |T] corresponding to one of the samples considered in 
the experiments in [3]. In this case, the particles have 
larger saturation magnetization (r = 247) but the dis- 
persion is more diluted (<f>Q = 4.64 x 10~ 6 ). We have 
performed Langcvin Dynamics simulations as well as sim- 
ulations employing our new methodology, as in Case 1. 
These simulations will be denoted as LD247 and CG247, 
respectively. A list of relevant simulation parameters for 
LD247 and CG247 simulations is given in Table |TT| The 
results obtained for (N(t)) are also given in Fig.TJT The 
results for the probability distribution p(s; t) at t = 5s are 
shown in Fig. [6] Again, we obtain a good agreement be- 
tween the predictions of both simulation methodologies, 
both in the average size of chains and in the probability 
distribution of chain sizes. 



As shown in Table [IT] both for the case with T = 40 
and r = 247, the computational cost of the CG simula- 
tion technique is much lower than the corresponding LD 
simulation. For example, we note that a production run 
of 6 seconds for the LD247 simulation requires about 4.5 
days of calculations, with the program running in parallel 
in a 8-Core AMD Opteron Magny Cours 6136 processor. 
In contrast, the CG247 simulation requires less than 4 
hours to simulate the same physical time using a single 
core of the same processor employed in the LD247 run. 
In addition, we can reach surprisingly long time scales 
in our CG247 simulation with a very low computational 
cost (see Table |TI| . In simulation CG247, we reach sim- 
ulated times up to 10 3 s in a 1 day calculation, a time 
scale two orders of magnitude larger than that accessible 
using Langevin Dynamics simulations. 

The CPU costs shown in Table [IJJ demonstrate that 
the new simulation technique allows us to perform simu- 
lations of the two systems considered here with an ex- 
tremely reduced computational effort as compared to 
Langevin Dynamics simulations. Moreover, it is also im- 
portant to notice that the required CPU time for the CG 
approach to simulate a certain time interval is reduced 
during a simulation, since the number of CG objects de- 
creases as the simulation advances. This effect is clearly 
shown in Fig. [7] for the CG40 simulation. We also show 
that the rate between the elapsed CPU time and the 
corresponding real time simulated depends linearly with 
the number of CG objects (see inset Fig. [7j. It should 
be noted that in the LD simulations the opposite effect 
was observed; the fact that the individual motion of the 
particles inside the chains are fully resolved makes LD 
simulations increasingly inefficient as time goes on. 
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FIG. 7: (Color online) Elapsed CPU time (open circles) and 
number of CG objects (solid circles) as a function of the real 
simulated time for the CG40 simulation. Inset: Rate be- 
tween the elapsed CPU time and the real simulated time as 
a function of the number of CG objects present in the CG40 
simulation. 



B. Further discussion on the approximations of the 

model 



1. Diffusion Model 

As it has been already mentioned, one of the two key 
ingredients of the CG approach is the diffusion model 
adopted to describe the motion of the coarse grain ob- 
jects. In order to check the possible influence of the 
model selected, we have also performed simulations with 
a different diffusion model proposed by Tirado et al. [28] 
in which they describe the translational motion of right 
circular cylinders also accounting for the so called end- 
effects. Following the same approach than in |23| . we 
have used the expressions: 



Dl _ 3 
~D[ ~ 2s 

D x 4s 

where 7" and 7 _L are the end-effect functions defined as: 

0.90 



= i-[ln(s)+j ± (s)}. 



(9) 
(10) 



7 N (s) = -0.21 



s 

0.18 0.24 



7 x (s) =0.84+ + 



(11) 
(12) 



We have computed the mean aggregate size (N) using 
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FIG. 8: (Color online) Comparison between the results ob- 
tained from the elongated rod approximation (correspond- 
ing to Eqs. Q and (JsJ and represented by solid circles) and 
the cylinder approximation with end-effects (corresponding 
to Eqs. (|9| and and represented by open circles) when 
computing the mean chain size of the aggregate (N) for the 
CG40 system. 



both diffusion models for the CG40 system and the re- 
sults obtained are plotted in Fig. [5J As it can be seen 
from these results, no significant differences are found in 
the average number of particles (N) obtained with both 
diffusion models. For this reason we can conclude that 
both models are suitable for the description of the diffu- 
sive motion of the chain-like aggregates in such systems 
and our selection of the elongated rod model (Eqs. Q 
and (5|) instead of Eqs. ([£])-( 12 1 for the simulations was 
based on its major simplicity. 



To this end, we have performed two additional simu- 
lations in which the s dependence of r a (s) is ignored. A 
first simulation (denoted as CG40-min) corresponds to a 
repetition of the simulation CG40 of the previous section 
(see Table [If]) but using r a = 1.46d for all chains. We also 
performed another simulation (denoted as CG40-max) in 
which we employed the value r a = 2.20c? for all chains. 
These values correspond to the minimum and maximum 
values of r a (s) employed in the original CG40 simulation 
(see Fig. [TJ. 

All these approaches give us different dynamics of the 
system as it is shown in Fig. [9] where the mean aggre- 
gate size (N) is plotted as a function of time together 
with the corresponding LD results. We observe that the 
CG40 simulation evolves from an initial behavior close to 
the CG-min simulation to a behavior closer to the CG- 
max simulation. Analogous calculations for the CG247 
system (not shown here) exhibit identical behavior. In 
consequence, is important to take into account the full 
r a (s) dependence in the simulations as described in our 
formulation of the model in Section II. 
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2. Effective interaction. Attraction radius 

As explained in detail in Section II, we have defined 
the aggregation regions for each CG object as the sur- 
rounding space in which the magnetic interaction energy 
between the CG object and a dummy single particle is 
equal or smaller than —ksT. As shown in Fig. [T] the 
attraction radius r a defining this region depends on the 
size of the considered aggregate and on the magnetic cou- 
pling parameter Y (see Eq. [6]). It is observed that for 
small chains the attraction radius increases with their 
size and tends to a constant value for larger chains (the 
addition of a new particle into the same aggregate does 
not significantly contribute to the interaction magnetic 
energy). Here, we would like to demonstrate the impor- 
tance of accounting for the s dependence of r a (s) in the 
simulations. 



FIG. 9: (Color online) Comparison of the calculated mean 
chain size (N) between two simplified versions of the coarse- 
grain model (see main text for details) and the full version. 
Crosses correspond to the CG40-min simulation, stars to 
the CG40-max simulation and solid circles correspond to the 
CG40 (full version) simulation. We also show the Langevin 
Dynamics results (LD40, open circles). 



C. An example of practical application: Chain 
growth and T2 measurements 

Our objective in this subsection is to illustrate the ap- 
plicability of the methodology developed here in situa- 
tions of interest for applications of superparamagnetic 
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TABLE III: Set of parameters used for numerical integration in the Coarse Grain simulations of the same colloids described 
in Case 2 in Table [I] but for different volume fractions <f>o. N p is the number of particles in the simulation, L z and L x — L y 
are the sizes of the simulation box (in units of particle diameters) in the directions parallel and perpendicular to the magnetic 
field, respectively. At is the time step and £/ is the total simulated time. As in Table |Tl} we also indicate the total amount of 
CPU time employed in the calculation. 
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FIG. 10: (Color online) Average number of colloids in a chain 
(TV) as a function of time for the simulations described in 
Table lUll 



colloids. As an example, let us consider the use of su- 
perparamagnetic colloids as contrasts agents in magnetic 
resonance imaging. An important issue in this applica- 
tion is the possibility of chain formation of colloids due 
to the strong magnetic fields applied in the experiments. 
The formation of chains of colloids in the sample increases 
the transversal relaxation time T 2 of protons, which is an 
undesired effect in practice. In a previous work, we have 
employed a preliminary version of our simulation code to 
analyze this possibility in MRI |9_. We have found that 
under conditions of interest for MRI, significant chain- 
ing occurs. We would like to discuss here the results of 
our simulations as well as compare our results with the 
experiments in a more direct way than the preliminary 
simulations presented in Ref. [9]. 

The system considered here is a dispersion of super- 
paramagnetic colloids in water with the physical proper- 
ties of Case 2 in Table U but now we have considered 4 
different values of the initial volume fraction of colloids. 



according to the experiments in Ref. [9] . The correspond- 
ing volume fractions arc given in Table |III| The simula- 
tions of these systems, performed with the methodology 
discussed in section II have been labeled as CG247-00, 
CG247-01, CG247-02 and CG247-03, respectively and all 
the technical details are given in Table |III| (note that the 
CG247-02 simulation in this table is identical to the sim- 
ulation CG247 of Table [n). 

The results for the average number of particles in a 



chain (N(t)) are given in Fig. 10 for the time scales rel- 
evant in the experiments of Ref. [9 . In all cases we ob- 
serve significant chain formation even in the case of the 
smallest concentration. Of course, the kinetics of chain 
formation is observed to slow down as the concentration 
of colloids decreases. 

The formation of chains has direct impact on the 
transversal relaxation rate 1/T 2 of water protons. Ini- 
tially (t — 0) , the relaxation rate of water protons 1 /T^ 
is determined by the presence of a random distribution 
of isolated (dispersed) colloids. As time goes on, chains 
form and modify the T 2 response of the surrounding wa- 
ter protons. Therefore, the experimentally measured T 2 
at a given instant t depends on the distribution of chain 
sizes at that time t. As proposed in Ref. [S], we can give a 
theoretical prediction of the relaxation rate 1/T 2 (f) from 
our simulation results by computing the following aver- 
age: 



1 



1 



T 2 (t) iV ? 



In Eq. (13 1, l/T^' is the relaxation rate of water pro- 



1 



T. 



(13) 



tons near a colloid forming part of a chain containing 
exactly s colloids and n s (t) is the number of chains of 
size s at time f, as defined in Section II. Our simulation 
results provide n s (t) whereas the calculation of l/T^ 
requires an additional study of the motion of water pro- 
tons near a chain containing s colloids. For the particles 
of the experiments (Case 2 in our Table [I]), the theoretical 

results for 1/T 2 (s) were given in Fig. 9 of Ref. [5]. These 
results can be well fitted to an analytical expression of 
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the form: 



T, 



T, 



(0)' 



(14) 



where a fit to the calculations in Ref. [5] gives a = 0.0415 
and b = 0.45. Now, making use of such a fit in Eq. [\ 
and the values of n s (i) obtained from simulations CG247 
00, CG247-01, CG247-02 and CG247-03, we can make 
a theoretical prediction for the relaxation rate 1/T<z(t). 
The results are compared in Fig. 11 with experimental 
results. The simulations show a remarkable agreement 
between theory and experiments for times corresponding 
to mean chain length (N) larger than 50 colloids. It 
should be noted that in the case of very long chains, 
the measurements are not reliable due to sedimentation 
effects 0. 
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FIG. 11: (Color online) Evolution of the relaxation rate I/T2 
of water protons in 4 dispersions containing different concen- 
trations of superparamagnetic colloids. Solid lines correspond 
to the predictions obtained from the simulations described in 
Table III and Eqs. (131 and (141. Symbols correspond to ex- 



perimental data extracted from Fig. 5a in Ref [9J. 



IV. CONCLUSIONS 

In this paper, we have presented a new on-the-fly 
coarse grain model to describe the chaining phenom- 
ena observed in dispersions of superparamagnetic colloids 
under strong external magnetic fields. We report sim- 
ulation results with the new methodology, which show 
good agreement with those obtained from more detailed 
Langevin Dynamics (LD) simulations. The great advan- 
tage of the new methodology presented here is its low 
computational cost in terms of CPU time. As a conse- 
quence, we are able to run longer simulations, reaching 



time scales not accessible in LD simulations. In order to 
illustrate the applicability of the code in experimentally 
relevant situations, we have considered the waiting time 
dependence of the relaxation rate I/T2 of water protons 
observed in magnetic resonance experiments of disper- 
sions of superparamagnetic colloids [9]. Experimental 
results corresponding to waiting times from 1 to 10 3 s 
were correctly predicted by our simulations. 

The model, in its present formulation, cannot be ap- 
plied to situations more complex than irreversible chain 
growth. However, it seems possible to expand the model 
to consider other situations of interest. A first generaliza- 
tion could involve the inclusion of lateral interactions be- 
tween the chains [29 , which are responsible for the forma- 
tion of thick chains, observed at volume fractions larger 
than those considered here [T7]. Optical microscopy ob- 
servations [TT] also show the formation of thick chains 
and bundles in magnetophoresis experiments (motion of 
magnetic particles under magnetic gradients). Hence, the 
inclusion of lateral interactions and deterministic motion 
of the CG objects will be needed in order to extend our 
model to study magnetophoresis. Another interesting ex- 
tension, which is now under way, is the inclusion of the 
possibility of chain breaking due to thermal fluctuations, 
a mechanism relevant at low values of the magnetic cou- 
pling parameter T. This extension of the model will allow 
us to study in depth the equilibrium state described in 
Refs. QI1 [30]. 

A final improvement to the model could be taking into 
account the full magnetic response M(H) of the particles 
in the simulation, in order to simulate situations in which 
the external magnetic fields are not strong enough to 
saturate the magnetic colloids. This is a typical situation 
in many published experimental studies of aggregation of 
magnetic colloids (see for example [231 124| ). which focus 
on the linear magnetic response regime of the colloids. 
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Appendix A: Selection of the integration timestep 

An important issue that one has to take into account 
when performing Brownian Dynamics simulations is the 
proper selection of the integration timestep. The typical 
diffusive displacement £ for a single colloid of diameter 
d and diffusion coefficient D after a timestep At can be 
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FIG. 12: (Color online) Influence of the integration timestep 
At in the simulations CG40 and CG247 (top and bottom, 
respectively) on the time evolution of (N). 



estimated by 



V6DM = d^QAt/r. 



(Al) 



where we have defined the characteristic diffusion time 
t = d 2 /D. In general, one selects a timestep At which re- 
sults in a displacement £ smaller than the relevant length 
scales of the problem (typical separations between par- 
ticles, range of interaction forces,...). In our model, the 
length scale of interactions is given by the radius of the 
attraction zones (see FigJTJ) . The typical diffusive dis- 
placement I corresponding to the selected At (Eq.(All) 



has to be smaller than the radius of the attraction zone 
r a of the CG objects. In this way, CG objects will cor- 
rectly explore the attraction zone of other surrounding 
CG objects. As it is shown in Figjl]:, the size r a of the 
attraction zone depends on the chain length s (the small- 
est value of r a corresponds to s = 1) and on the coupling 
parameter T. The dependence on T is strong, so one has 
to take into account this fact in selecting At. 

In order to check the effect of At in the results of our 
simulations at T — 40 and V = 247, we have repeated 
the CG40 and CG247 simulations with three different 
timesteps: At = l.OOr, O.IOt and O.Olr. These timesteps 
At correspond to typical displ acem ents of £ ~ 2.4c?, 0.8c? 
and 0.24c? respectively (see Eq.( Al )). The results of these 
simulations for the average number of particles in a chain 
are shown in Fig 12 As it can be observed, the effects 
of the selection of the timestep are critical for the CG40 
system and irrelevant for the CG247 system. In order to 
understand the effect of these different At , one has to 
compare the I obtained for each At with the values of 
r a (s = 1) calculated for r = 40 and T = 247 (see Fig- 
ure [TJ:). In the CG40 simulation, the smallest radius of 
a CG attraction zone is r a (s = 1) = 1.46c? (see Figure 
[lp, case r = 40), so the attraction sphere has a diam- 
eter similar to the displacement £ = 2.4c? obtained for 
At = l.Or. In this case, colloids cannot explore properly 
the attraction zones and many chain formation events are 
lost during the simulation run. The other two At give al- 
most identical results since in both cases I is smaller than 
T" a {s — 1), and attraction zones are explored properly. In 
the CG247 simulations, we observe almost identical re- 
sults for the three selected At. In this case, we have 
r a (s = 1) = 2.89c? (see Figure [TJ;, case T = 247), which 
is larger than the corresponding typical displacements I. 
As a general rule, if the timestep selected is too large, 
the chains cannot properly explore the binding sites of 
the surrounding chains and the chaining process is not 
correctly simulated. 
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